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Abstract. We analyze 2.8-yr data of 1-100 GeV photons for clusters of galaxies, collected 
with the Large Area Telescope onboard the Fermi satellite. By analyzing 49 nearby massive 
clusters located at high Galactic latitudes, we find no excess gamma-ray emission towards 
directions of the galaxy clusters. Using flux upper limits, we show that the Fornax cluster 
provides the most stringent constraints on the dark matter annihilation cross section. Stack- 
ing a large sample of nearby clusters does not help improve the limit for most dark matter 
models. This suggests that a detailed modeling of the Fornax cluster is important for setting 
robust limits on the dark matter annihilation cross section based on clusters. We therefore 
perform the detailed mass modeling and predict the expected dark matter annihilation sig- 
nals from the Fornax cluster, by taking into account effects of dark matter contraction and 
substructures. By modeling the mass distribution of baryons (stars and gas) around a cen- 
tral bright elliptical galaxy, NGC 1399, and using a modified contraction model motivated by 
numerical simulations, we show that the dark matter contraction boosts the annihilation sig- 
natures by a factor of 4. For dark matter masses around 10 GeV, the upper limit obtained on 
the annihilation cross section times relative velocity is (av) < (2-3) x 10~ 25 cm 3 s _1 , which 
is within a factor of 10 from the value required to explain the dark matter relic density. 
This effect is more robust than the annihilation boost due to substructure, and it is more 
important unless the mass of the smallest subhalos is much smaller than that of the Sun. 
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1 Introduction 

Revealing the nature of dark matter is one of the important and challenging goals of modern 
physics, astrophysics, and cosmology. One of the avenues involves searching for annihilation 
signatures from dark matter halos in a form of gamma rays, neutrinos, and charged particles. 
Such a possibility is theoretically well motivated especially if dark matter is made of weakly in- 
teracting massive particles (WIMPs) such as the super symmetric neutralino [1]. In addition, 
in order to thermally produce dark matter at the right abundance in the present Universe, 
WIMPs have to have some certain annihilation cross section (times relative velocity), most 
preferably, (av) = 3 x 10 -26 cm 3 s _1 [1-3]. The recently launched Fermi Gamma-Ray Space 
Telescope and ground-based atmospheric Cerenkov telescopes such as HESS, MAGIC, and 
VERITAS have started to provide interesting constraints on the annihilation cross section of 
dark matter particles from various astrophysical objects, ranging from dwarf galaxies to the 
isotropic diffuse gamma-ray background, i.e., the emission coming from the entire Universe. 

Clusters of galaxies are the largest concentrations of dark matter pulled together by 
gravity, and therefore, an attractive astrophysical object to look for the annihilation signals. 
The Fermi Large Area Telescope (LAT) collaboration analyzed their 11-month data to place 
constraints on the dark matter annihilation cross section [4] and decay [5] using observations 
of several nearby galaxy clusters. The obtained upper limits have already started to exclude 
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some of the parameter space of WIMPs, and these cluster-based constraints are quite compa- 
rable to those obtained from dwarf galaxies [6-8] and the diffuse gamma-ray background [9] . 

In the previous analyses as well as theoretical studies on dark matter annihilation in 
galaxy clusters, it is commonly assumed that the dark matter density profile is smooth and 
characterized by the Navarro- Frenk- White (NFW) form [10]. However, it has been pointed 
out that the effects of subhalos might significantly boost the dark matter annihilation signal 
for massive clusters [11-13], if the smallest subhalo masses are of around Earth size as implied 
by arguments of the kinetic decoupling of WIMP particles in the early Universe [14-18]. In 
the case of the smooth NFW-like halos, the current upper limits on the cross section is about 
~30 times larger than the canonical value ({o~v) = 3 x 10~ 26 cm 3 s" 1 ) for the low-mass dark 
matter particles around ~10 GeV. In the case of Earth-mass substructures, on the other 
hand, the canonical cross section has already been excluded for the low-mass WIMPs [4, 19], 
although there is a large uncertainty on the subhalo models, especially since it involves 
significant extrapolation of the subhalo mass functions down to small scales. 

There is another effect that must be taken into account but has received little attention 
to date — the effect of dark matter contraction due to baryonic infall. Unlike dark matter, 
baryons can lose energy via radiative cooling and they fall towards the center of the halo, 
and as the result, modify gravitational potential in the central region. Dark matter is then 
dragged towards the center because of the deepened potential, modifying the density profile 
of dark matter [20-23]. Its implications for dark matter annihilation have been discussed in 
the context of galaxy-size halos [24, 25], but not in the context of cluster-size halos. Galaxy 
clusters contain a large amount of baryonic gas and stars, and in fact, the stars dominate the 
gravitational potential at central regions. Although there is no firm observational evidence 
for or against the halo contraction in massive clusters at the moment, the effect of the dark 
matter contraction is likely important for cluster-size halos. 

In this paper, we analyze Fermi-LAT 2.8-yr data in order to constrain the dark matter 
annihilation signatures from galaxy clusters. As we shall show in this work, among several 
clusters analyzed thus far, the most stringent constraint on the annihilation cross section 
is obtained with Fornax, a nearby, well measured galaxy cluster located at 20 Mpc away 
from Earth. We therefore focus on interpreting observations of the Fornax cluster, paying 
particular attention to modeling of dark matter contraction due to baryon dissipation as well 
as substructure boost, and address the relative importance of these two effects. Apart from 
its proximity and its large mass (~lO 14 /t _1 M ), there are several reasons why the Fornax 
cluster is ideal for detailed modeling. First, the Fornax cluster does not host any bright 
active galactic nuclei, unlike the Virgo and Perseus clusters. Secondly, it has regular thermal 
gas profiles and a spherical central massive elliptical galaxy, NGC 1399, and this feature 
makes the calculation of contraction based on the assumption of spherical symmetry quite 
reasonable. In order to infer the distribution of baryons, we use high-resolution observations 
of surface brightness profiles with Hubble Space Telescope (HST) for stars and with ROSAT 
X-ray Telescope for thermal gas. Using the well-measured stellar and gas profiles, we then 
compute the dark matter density profile taking into account the dark matter contraction. 

This paper is organized as follows. In Sec. 2, we present the formulation to compute 
intensities of gamma rays from dark matter annihilation, and models of dark matter density 
profiles, halo contraction due to baryon dissipation, and substructures from numerical sim- 
ulations. In Sec. 3, we discuss the results of analysis for a large sample of galaxy clusters, 
including the upper limits on individual clusters and stacking analysis, highlighting that a 
detailed mass modeling of the Fornax cluster is indeed more important than the statistical 
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constraint. Procedures for the data analysis of the Fermi-LAT data are summarized in Ap- 
pendix A. In Sec. 4, we show results for the Fornax cluster, in which we present a detailed 
mass modeling of dark matter, stellar, and gaseous components in Fornax and discuss the 
relative importance of halo contraction and subhalo boost effects in modeling the annihila- 
tion signal from clusters. We finally conclude the present paper in Sec. 5. Throughout the 
paper, we assume the flat cold dark matter cosmology with cosmological constant (ACDM), 
and adopt U m = 0.277, A = 0.723, H Q = lOO/ikms" 1 Mpc" 1 , and h = 0.7. 



2 Gamma-ray intensity and dark matter density profiles 



In order to compute the gamma-ray intensity from dark matter annihilation, one needs to 
specify the dark matter annihilation mechanisms and the dark matter distribution in clusters. 
Here, after briefly introducing formalism to compute the intensity in Sec. 2.1, we discuss 
the density profiles obtained with numerical N-body simulations in the literature, where 
gravitational structure formation and evolution of only dark matter particles are followed. 
We introduce the smoothly distributed large-scale halo component (Sec. 2.2), effects of halo 
contraction due to baryon dissipation (Sec. 2.3), and small-scale substructure (Sec. 2.4). 

2.1 Gamma-ray intensity from dark matter annihilation 

The gamma-ray intensity, here defined as the number of photons received per unit area, time, 
solid angle, and energy, i.e., i~ = dN^ / dAdtdO.dE , from annihilation of dark matter particles 
X, in a halo is given by 



(av) 



dN, 



7,ann 



8vr(l 



dE' 



E'=(l+z)E 



ds p 2 [r(s,6)\, 



(2.1) 



where is the angle from the halo center, z is the redshift of the halo, m x is dark matter mass, 
and d-ZV^ann/ 'dE is the gamma-ray spectrum per annihilation. The intensity is proportional to 
the line-of-sight (s) integral of dark matter density squared, p^. Here, the density is assumed 

to be spherically symmetric, and the radius r is related to s and 9 through r = \J s 2 + d? A 6 2 , 
where d& is the angular diameter distance to the halo. 

The gamma-ray spectrum depends on the final state of annihilation. In this study, we 
adopt three different annihilation channels, XX ~> bb, W + W~, and t + t~~. Figure 1 shows 
E 2 dN^ ann /dE for various dark matter masses between 10 GeV and 1 TeV. The annihilation 
spectra for neutralino WIMPs are typically given by some combination of these three spectra. 
Note that there is no annihilation happening if the final state is W + W~ and dark matter 
mass is smaller than mass (80.4 GeV). 

2.2 Dark matter density profile: smooth component 

For the smooth component of the dark matter density profile, p(r), we consider two models: 
one is the NFW profile [10], and the other is the Einasto profile with a central core [26, 27]. 
The NFW profile is given by 



p(r) 



(2.2) 



(r/r s )(r/r s + l) 2 > 

where p s and r s are the scale density and scale radius, respectively. The virial mass M v \ r 
and redshift z are given inputs, and. the virial radius r v j r is obtained by solving M v - ir = 

; and 
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c(z)p c (z)/3, where p c (z) is the critical density, A v j r (z) = 18-7T + 82d — 39d 
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Figure 1. Gamma-ray spectrum per dark matter annihilation, E 2 dN 7tllnn /dE for (a) bb, (b) W + W~ , 
(c) t + t~ annihilation channels. Assumed dark matter masses arc m x = 10, 25, 50, 100, 250, 500, 
and 1000 GeV for the panels (a) and (c), and the same unless m x < mw = 80.4 GeV for the panel 
(b). 



d = Q m (l + z) 3 /[Q m (l J t- z) 3 + £1^] — 1 [28]. The scale radius r s is then denned as r s = r v ; r /c v i r , 
where c v ; r is the concentration parameter. We adopt the result of the mass-concentration 
relations from Ref. [29] (see also Refs. [30-33]): 



C v i r (iWvi r , z) 



7.85 



M v 



(1 + z) - 71 \2x lO 12 ^ 1 !/, 







-0.081 



(2.3) 



By taking volume integral of p(r) up to the virial radius r v j r and then equating it to M v i r , 
we obtain 



Airr 3 
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The Einasto profile, on the other hand, has the following functional form: 




l/n 



p(r) = p e exp < -d, 



- 1 



(2.5) 



where n and r e are the free parameters, and d n = 3n — 1/3 + 0.0079/n. The normalization 
p e is obtained such that the total mass within r v j r agrees with the measured M v i r : 



where T(a, x) is the lower incomplete gamma function. 

Note that we also use other definitions of mass, including Ma with A = 200 or 500, 
defined as the enclosed mass within a sphere with radius r a , within which the average density 
is A times the critical density of the universe. Conversion between different mass definitions is 
straightforward if one assumes density profile such as NFW as well as the concentration-mass 
relation (e.g., [34]). 

2.3 Dark matter contraction due to baryon dissipation 

The dark matter density profile is modified as a result of baryon dissipation. A so-called 
dark matter contraction model is based on the conservation of the angular momentum: 



where M(r) is the mass enclosed within the radius r, subscripts dm and b are for dark 
matter and baryons, respectively, and subscripts i and / represent quantities before and 
after the contraction, respectively. The masses inferred from, e.g., gravitational lensing or X- 
ray observations are the total mass, and therefore one must assume underlying distributions 
of both baryons and dark matter in order to infer -Md m ,i(r). Here, when we discuss the 
initial density profile (before contraction), we simply assume that dark matter and baryons 
trace each other exactly; namely, Pdm,i(r) = (1 - fb)Pi(r) and M dmii (r) = (1 - f b )Mi{r), 
where /& is the baryon fraction at virial radius. For the final profiles of baryons, on the other 
hand, we use the observed distribution of stars and gas (as discussed in the later section), 
giving the baryon mass distribution to be M b j{r) = M±(r) + M gas (r). Finally, f is the orbit- 
averaged radius for a particle found at r in a simulation [22]. Equation (2.7) is solved for r/, 
a contracted radius, as a function of the initial radius rj. The dark matter mass profile after 
the contraction is then obtained by M& m j(rf) = M& m ^(ri{r f)) . 

If particle orbits can be approximated to be circular and if the contraction happens 
adiabatically, then we simply have f = r. This is a so-called standard adiabatic contraction 
model originally proposed in Ref. [20]. Reference [22], based on hydro dynamical simula- 
tions, presented the modified formula (2.7) and pointed out that corrections to the adia- 
batic contraction model can be accommodated by modifying a relation between f and r 
as f/r v i r = 0.85(r/r v i r ) ' 8 . This was further investigated by various groups both theoreti- 
cally [25, 35-41] and observationally [42-45]. Recently, Ref. [23] systematically studied this 
effect, and found that the relation 




(2.6) 



Miir^ri = [M dmji (fi) + M bJ (r f )]r f 



(2.7) 




(2.8) 



- 5 - 



well represents the simulation results, where ro = 0.03r v i r , and w ranges from 0.5 to 1.3, 
if Aq is fixed to 1.6. Here we adopt w = 0.8 as a canonical value for this parameter, 
but also consider the cases of w = 0.6 and 1. In addition, we adopt Aq = 1 and w = 
1, the standard adiabatic contraction, as an extreme case scenario. Equation (2.7) has 
been calibrated only down to r ~ 10~ 3 r v i r because of the resolution limit [23], but we shall 
extrapolate this relation further to very small radii. Although we assume angular momentum 
conservation, there is no guarantee that Eq. (2.8) remains valid below the resolution limit 
of the simulations. We therefore adopt yet another relation as a conservative scenario of 
the contraction: f = maxLAoro^/ro)™, 10 _2 r v i r ], which does not allow f to reach below 
10 _2 r v ; r , the current resolution limit of simulations. Note that in this model the magnitude 
of contraction is smaller for the larger f. We discuss the application of this model to the 
Fornax cluster in Sec. 4. 

2.4 Effects of substructure 

If galaxy clusters consist of numerous subhalos as implied by numerical simulations, then the 
gamma-ray radiation may be dominated by the contribution from subhalos, even though they 
make up a small fraction of the total cluster mass. Here, we consider the results of recent 
numerical simulations of cluster-size halos that contain subhalos of M su b > 5 x 10 7 Mq, and 
its extrapolation down to Earth mass, M m \ n = 10 -6 M Q , or even smaller mass 10 _12 Mq [12]. 
(This is because the minimum mass of subhalos for neutralino dark matter models ranges 
widely, M Ui \ n = 10~ 12 -10~ 4 Af Q [18].) According to this model, the surface brightness profile 
contributed by subhalos are given by 

T , m _ 6(M 200 )F 7 , host (ff) 1 

7 ' subi ' J_ 7rln(17) 2 + (W4) 2 ' 1 ' 

/ yr \ 0.39 

6(M 200 ) = 1.6 x 10~ 3 f-^-j , (2.10) 

for 6 < 6200 and for M m ; n = 1O _6 M0. Here #200 = ^ooAUi the flux F 7 is obtained by 
integrating I 7 [Eq. (2.1)] over the solid angle, b = i^sub/^host is t ne boost factor, and the 
subscripts "sub" and "host" indicate subhalos and the host halo, respectively. Note that the 
boost factor decreases significantly for the larger values of minimum subhalo mass M m i n , and 
it scales roughly as b oc M~® n 2 [12]. 



3 Constraints from a large statistical sample of clusters 

We discuss constraints on annihilation cross sections from 49 nearby rich galaxy clusters. 
After computing the intensity I^(E,9) with Eq. (2.1) for these clusters, we analyze the 
Fermi-LAT data for 2.8 years to place constraints on (crv). Unlike the previous study where 
clusters were assumed to be point-like sources [4], here we take the extension of the cluster 
emission fully into account. Analysis details can be found in the Appendix A. In this section, 
we consider only the NFW smooth halo model (Sec. 2.2) for a density profile, and other 
models are discussed in the following section. 

3.1 Upper limits on individual clusters 

We choose 23 galaxy clusters that are the brightest in X-rays from Ref. [46]. These are the 
same clusters as those analyzed in Ref. [47] for general gamma-ray analyses of galaxy clusters, 
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Table 1. Parameters of galaxy clusters analyzed for dark matter annihilation. Upper sample is from 
Chandra Cluster Cosmology Project [48] and lower one is for the X-ray-selected bright clusters [46]. 
A2199 and A3571 arc found in the both samples. 
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where the sample clusters are chosen such that they yield the largest mass divided by distance 
squared, which is a reliable measure of expected annihilation fluxes. In addition, we also 
analyze 34 massive low-redshift clusters from the Chandra Cluster Cosmology Project [48]. 
We include this sample in order to investigate the level of cosmic-ray pressure in clusters 
and their implications for cluster-based cosmological tests [49] , which we will discuss in more 
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Figure 2. The cross section upper limits from Fornax cluster (smooth NFW with no contraction) 
for the bb, W + W~ , and t + t~ annihilation channels. 



detail in our subsequent work. We remove clusters located at low Galactic latitudes, |6| < 20°, 
in the analysis because of the strong emission from the Galactic plane. Consequently, we 
have 49 clusters to analyze. Two of them (A2199 and A3571) are found in both samples. 
Table 1 summarizes properties of these clusters. 

In the table we also show the flux upper limits for each cluster assuming E~ 2 spectrum 
and that the clusters are point-like. The seventh column shows the integrated flux upper 
limits in 0.1-100 GeV, F Um (> 0.1 GeV), and they are typically (1-5) x 10 -9 cm^s -1 . The 
last column of the table shows the ratio of expected gamma-ray flux from dark matter 
annihilation in the case of smooth NFW halo model to the flux upper limits, with respect to 
the value for the Coma cluster. These ratios, therefore, indicate the promise of each cluster 
as a target for dark matter annihilation. It is clear from this table that the Fornax cluster is 
by far the most promising target among 49 clusters considered in this work. 

In Fig. 2, we show cross section upper limits as a function of WIMP mass for the three 
different annihilation channels obtained from observations of the Fornax cluster. The channel 
into gauge bosons will give similar yield to bb for all the masses, but limits for r-lepton pairs 
are significantly weaker especially at larger masses, because of hard spectrum (Fig. 1). We 
note that these limits are consistent with those of Refs. [4, 19] after correcting for differences 
in exposures. 

3.2 Stacking analysis 

One may wonder if stacking a large statistical sample of galaxy clusters improves upper limits 
on the annihilation cross section from those obtained individually. Here, we stack a subset 
of 49 galaxy clusters discussed above, which is a significant improvement from the previous 
studies based on 6 or 8 clusters [4, 19]. The procedure for obtaining the stacking constraint 
on the cross section is similar to the individual cluster analysis (see Appendix A for more 
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Figure 3. Upper limits on annihilation cross section (av) for the bb, W + W~ , and t + t~ annihilation 
channels, in the case of smooth NFW host halo model (with no contraction). Solid and dashed curves 
are the limits due to staked clusters and Fornax, respectively. 



details). In order not to bias the limits too high, we remove clusters with signal-to- noise ratio 
larger than 3cr compared with the (fixed) background model for any dark matter masses and 
annihilation channels. 1 As the result, we have 38 clusters used in the stacking analysis. 

In Fig. 3, we compare the cross section upper limits from the stacking analysis and 
the constraint from the Fornax cluster alone for the three different annihilation channels. 
We find that the stacked limits of the cross section degrade toward lower dark matter mass, 
while they improve by about a factor of 2 at multi-TeV mass regime for the t + t~ channel. 



1 Note that this does not mean that some clusters are already detected with > 3a, because there are 
uncertainties in our procedure of fixing background. 
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This is in part because we still have weak (less than 3<r) excess emission from some galaxy 
clusters for some annihilation models. This shows that the Fornax cluster dominates the 
constraint on the annihilation cross section, and stacking a large sample of clusters does not 
help improve the limits. This is in contrast to the case of dwarf satellite galaxies, where the 
stacking analysis of 10 satellite galaxies can improve the limits by a factor of a few for a wide 
range of dark matter masses [7, 8]. In the next section, we shall discuss that treating the 
Fornax cluster carefully is indeed more important for setting an improved constraint on the 
annihilation cross section based on galaxy clusters. 

4 Constraints from the Fornax cluster 

As shown in Table 1, the Fornax cluster is located at (I, b) = (236.72°, —53.64°) in the Galactic 
coordinate. Its redshift is z = 0.0046, and the corresponding angular diameter distance is 
(1a = 13.7 h~ l Mpc. Both optical and X-ray observations show that the brightness profiles of 
Fornax are fairly regular, and therefore, modeling dark matter distribution assuming spherical 
symmetry is well justified. 

4.1 Distribution of dark matter, stars, and gas 

The X-ray brightness has been used to infer the mass to be M500 = 6.4 x 10 13 h~ l M Q [46, 50] 
(see also Refs. [51, 52]), and the corresponding virial mass and concentration are M v j r = 

1.2 x 10 14 h~ l Mq and c v ; r = 5.6, respectively. For the Einasto profile, Eq. (2.5), we adopt 
n = 3.87 and r e = 1247 kpc, following the cluster simulation of similar size (model C09 of 



We shall obtain M*(r) and M gas (r) as they are necessary ingredients for calculating the 
final dark matter profile after contraction using Eq. (2.7). For the stellar mass distribution, 
we use the HST observation of the bright elliptical galaxy, NGC 1399, located at the center of 
the Fornax cluster. The surface brightness profile of stars are well represented by the "Nuker 



which behaves as 9 7 (9 ") at small (large) 9 with a break around 9^. We here consider the 
results of HST WFPC2 observations using the F606W filter [54, 55]: a = 1.58, (3 = 1.63, 
7 = 0.09, 9\) = 3.17", and 1^ = 17.23 mag arcsec -2 . This photometry result is valid down to a 
very small radii, ~1 pc (i.e., ~0.0l"). The same measurements also constrain the ellipticity of 
the galaxy, which is found very small [54] . Given that the assumption of spherical symmetry 
is well justified, we obtain the luminosity density by de-projection [56]: 



In order to obtain the density profile, we finally multiply j(r) by the mass-to-light ratio, 
P*( r ) = {M/L)j(r). According to the stellar kinematics study, it is believed that the gravi- 
tational potential in the central regions is dominated by stellar component, and the mass-to- 
light ratio has been inferred to be about 10 times that of the Sun [57, 58], which is typical 
for old stellar population. Figure 4 shows the density profile of the stellar component, p*, 
obtained in this way, and it dominates the mass within ~10kpc. We note that the power-law 



Ref. [26]). 



law" [53]: 




(4.1) 




(4.2) 
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Figure 4. Density (left) and mass (right) profiles of dark matter (before contraction in short-dashed 
line; after contraction in solid line), stars (long-dashed), and baryonic gas (dotted). For dark matter 
before and after contraction, both the NFW (thick) and Einasto (thin) profiles are shown. The 
canonical parameters (Aq = 1.6 and w = 0.8) are adopted for the contraction model. 



extrapolation of Eq. (4.1) works reasonably well even at large radii, but underestimates the 
brightness (and therefore gives conservative results) beyond ~20kpc [59]. 

The thermal gas component can be inferred from X-ray data. Here we use the ROSAT 
observations of the Fornax cluster [60, 61]. In particular, Ref. [61] identified three different 
components and fitted each with a bi-dimensional /3-profile. The thermal gas density and 
mass profile are shown in Figure 4. They are subdominant compared to the stellar density and 
mass within ~100kpc, and hence do not play important role in the dark matter contraction 
directly. However, since it is the dominant baryon component around the virial radius, it 
affects the dark matter contraction calculation through in an indirect way. 

4.2 Effects of dark matter contraction 

In Fig. 4, we show both initial and final dark matter density profiles for the Fornax cluster. 
We consider the NFW and Einasto models as the initial dark matter profile for the Fornax 
cluster and adopt the canonical contraction model (Aq = 1-6 and w = 0.8). The initial 
dark matter profiles are in good agreement with each other at radii larger than ~10 kpc, 
beyond which the total mass has been measured. Although the dark matter density is still 
subdominant at the smaller radii compared to the baryon density dominated by stars, it is 
enhanced significantly because of the dark matter contraction due to baryonic infall. This is 
the case even if the original profile has an inner core as the Einasto profile. 

Once the density profile is obtained, we compute the gamma-ray intensity using Eq. (2.1). 
Figure 5 shows the intensity integrated over energy, 1^(0), for both the NFW and Einasto pro- 
files as initial density models. Here, we fix particle-physics parameters to be Nj i3iXm (av) /m 2 = 
10~ 26 cm 3 s" 1 GeV~ 2 . In order to see what range of 6 contributes to the total signals the 
most, we multiply 1^(6) by 6 2 , such that the area below each curve gives the gamma-ray 
flux. One notices that there are double-peak structures. One peak at higher 9 corresponds 
to the scale radius r s of the NFW or Einasto profiles. The other peak, on the other hand, 
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Figure 5. Gamma-ray intensity as a function of angle from the cluster center before and af- 
ter the dark matter contraction, for both the NFW and Einasto profiles. Note that the instru- 
mental responses have not been included yet. The particle-physics parameters are fixed to be 
N 7 ^nn{crv)/m 2 x = lCT 26 cm 3 s- 1 GcV- 2 . 



corresponds to the characteristic scale of the dark matter contraction. Since the typical an- 
gular resolutions of the current generation of gamma-ray telescopes are on the order of 0.1-1° 
(depending on energies of gamma rays), the peaks at lower 9 will simply make a central pixel 
brighter but do not make additional resolvable features in the maps. 

The left panel of Fig. 6 shows a density enhancement factor after dark matter contrac- 
tion, Pdmj I Pdm,i, for various dark matter contraction models. Here we assume the NFW 
model as an initial density profile. We compare results for three different values of w = 0.6, 
0.8, and 1 by fixing Aq = 1.6 in Eq. (2.8). The larger the parameter w is, the more efficient 
the contraction gets, since f is smaller. We also show results of two extreme scenarios. One 
is the standard adiabatic contraction model (Ao = 1 and w = 1), and the other is the broken 
power-law model where f follows Eq. (2.8) with w = 0.6 at large radii but becomes constant 
(f = 10~ 2 r v ; r ) for small radii [see discussions after Eq. (2.8)]; they are labeled as "Adiabatic" 
and "Break," respectively. The true model is likely bracketed by these extreme cases. 

The right panel of Fig. 6 shows the gamma-ray intensities for the various contraction 
models considered above. As shown in the left panel of Fig. 6, the flux enhancements are 
bracketed by the two extreme models. Although the peak moves towards larger 9 for smaller 
f, one cannot resolve such peak yet with the current gamma-ray telescopes. The inner peaks 
are around 10 _5 -10 -3 deg, which corresponds to 3-300 pc in the physical size. By integrating 
over 9 out to large values, we find that the gamma-ray flux is boosted by a factor of 2-7 by 
the dark matter contraction effect. 

We show, in Fig. 7, the upper limits (corresponding to 95% credible interval) on the 
annihilation cross section, {crv), as a function of dark matter mass for the initial NFW and 
Einasto profiles, and for various annihilation channels. We compare no contraction model 
with the canonical contraction model (Aq = 1.6 and w = 0.8). The contraction model 
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Figure 6. Density enhancement of dark matter, p x (f)/ pNFw(f) (left) and gamma-ray intensity as a 
function of angle from the cluster center (right), due to the contraction for the initial NFW model. 
Solid curves are for the power-law models off [Eq. (2.8)], with w = 0.6 (top at low radii), 0.8 (middle), 
and 1 (bottom). Dashed line is for the standard adiabatic contraction model (Aq = l,w = 1), and 
dotted line is for A = 1.6, w = 0.6 but with a power-law break in Eq. (2.8) such that the f gets 
constant at 10 _2 r v i r for small radii. 



improves the limits at m x = 100 GeV by a factor of 4.1 for the bb and W + W~ channels, and 
3.8 for the t + t~ channel in the case of the NFW initial profile, and this is almost independent 
of mass. In the case of the Einasto initial profile, on the other hand, the improvement factors 
at 100 GeV are 2.4 for bb and W + W~, and 2.2 for r + r~. We note that the limits for no- 
contraction models of the NFW and Einasto profiles differ by only ~20% (the former is 
better), which is negligible compared to the differences between models with and without 
contraction. 

We compare various contraction models for the NFW initial profile and the bb annihila- 
tion channel in Fig. 8. As expected, the improvement factor is bracketed by our two extreme 
models, between 1.9 and 6.3 for the mass of 100 GeV and the bb channel. It is worth point- 
ing out that even the very conservative model of the contraction with the power-law break 
at small radii in Eq. (2.8) gives a factor of 2 enhancement compared to the no-contraction 
model, highlighting the importance of dark matter contraction in interpretation of the dark 
matter annihilation signals from galaxy clusters. This effect will make clusters brighter in a 
manner independent of the boost due to subhalos, which is the effect of much interest in the 
past several years and revisited below. 

4.3 Effects of substructure boost and minimum mass 

If the gamma rays are dominated by annihilation in subhalos and the subhalo mass function 
extends down to Earth mass (i.e., M min = 1O" 6 M ), then the cross section upper limits 
are significantly improved. Figure 9 shows the intensity 1^(9) integrated over energy, as a 
function of angle 9, for the Fornax cluster. Here, we fix the particle-physics parameter to be 
Nj >ami (av} / = 3 x 10~ 30 cm 3 s" 1 GeV" 2 . We show the contributions from both the smooth 
NFW host halo and substructures. For the latter, we adopt two minimum subhalo masses: 
5x 10 7 M Q corresponding to the current resolution limit [12] and 10 _6 Mq corresponding to the 
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Figure 7. Upper limits (corresponding to 95% credible interval) on annihilation cross section (av) as 
a function of dark matter mass. The limits for the initial profiles (no-contraction model) are compared 
with those for the canonical contraction model with Aq = 1.6 and w = 0.8. Thick and thin curves 
correspond to the NFW and Einasto profiles, respectively. 



typical smallest subhalo mass for the neutralino dark matter. One can see that the subhalo 
contribution could exceed that of the host halo if the subhalo mass function extends down to 
the Earth-mass scale. However, the relative importance of the host vs. subhalo contribution 
is highly sensitive to the cutoff mass scale. For example, if it is around the current resolution 
limit of simulations, then the host halo contribution becomes more dominant. 

Next, we investigate the dependence of the upper limit of dark matter annihilation 
cross section from Fornax on the assumed minimum subhalo mass M m ; n . Results are shown 
in Fig. 10 for m x = 100 GeV and the three different annihilation channels. Here we adopt 
the NFW profile with no contraction for the smooth component. Assuming the b oc M~j^ 2 
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Figure 8. Cross section upper limits (crv) for various contraction models for the NFW initial profile 
and the bb annihilation channel. The line types are the same as those in Fig. 6. 




Figure 9. Intensity 1(0) as a function of angle 9 (and radius in units of r v ; r ) for the Fornax cluster. 
Dark matter parameters are fixed to be N 7 , ann (<rv) /m^ = 3 x 10 cm s . Solid curve is the 
contribution from smooth NFW component, and dashed curves are those from subhalos assuming 
minimum subhalo mass of M m - m = 5 x 10 7 M Q (lower) and M m i n = 1O -6 M0 (upper). Total intensities 
are shown as dot-dashed curves. 
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Figure 10. Upper limits on annihilation cross section (crv), obtained with the Fornax cluster, as a 
function of minimum subhalo mass, M m j n , for the three different annihilation channels and m x — 
100 GeV. For the smooth component, the NFW profile with no contraction is adopted. 



dependence, we find that the upper limits on (av) scales as M^f n , when M m - m is small 
and the emission is dominated by annihilation in subhalos. This is realized for M m ; n < 
10 2 Mq. For larger M m j n than this value, the emission is dominated by the smooth NFW 
host-halo component and the limits asymptotically approach the values evaluated assuming 
no substructure. The canonical cross section of 3 x 10~ 26 cm 3 s" 1 is excluded, only if M m [ n 
is smaller than the Earth mass. We once again stress that the estimate of substructure 
boost factor (by Earth-size dark matter subhalos) requires extrapolation of the subhalo mass 
function by more than 13 orders of magnitude, given that the state-of-the-art numerical 
simulations can resolve halos of only M ~ 5 x 10 7 Mq level [12]. 

Comparing the relative importance of substructure boost to the effect of dark matter 
contraction, we conclude that the substructure boost could exceed the contraction boost, 
only if the minimum subhalo mass is considerably smaller than one solar mass. Otherwise, 
the substructure boost is at most a factor of a few, and the effect is generally smaller than 
the effects of contraction that ranges between a factor of 2 and 6. While the substructure 
boost remains highly uncertain, it is critical to take into account the contraction by baryon 
dissipation when studying the dark matter annihilation signals from clusters. 

5 Conclusions 

We investigated dark matter annihilation in galaxy clusters and constraints on annihilation 
cross sections with gamma-ray data from Fermi-LAT. Our main results are summarized 
below. 

By analyzing 49 nearby massive galaxy clusters located at high Galactic latitudes, we 
show that the Fornax cluster provides the most stringent constraints on the dark matter 
annihilation cross section. In addition, we show that the stacking analysis does not help 
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improve the cross section upper limits much. Our results suggest that detailed mass modeling 
of the Fornax cluster is more important for improving the dark matter annihilation cross 
section constraints from clusters than improving the cluster sample size. 

We therefore performed a detailed mass modeling and predicted expected dark matter 
annihilation signals of the Fornax cluster, by taking into account effects of dark matter con- 
traction and substructures. While the latter has been considered in the literature, the effect 
of dark matter contraction by baryon dissipation has not been considered in interpretation 
of the dark matter annihilation signal in cluster-size halos to date. By modeling the mass 
distribution of baryons (stars and gas) around a central bright elliptical galaxy, NGC 1399, 
and using a modified contraction model motivated by numerical simulations, we show that 
the dark matter contraction boosts annihilation signatures by a factor of ~4. Within a con- 
servative range of contraction parameters, we also show that the flux boost will fall between 
a factor of 2 and 7. 

We analyzed the Fermi-LAT data for 2.8 years around the Fornax cluster. After tak- 
ing all the sources and diffuse backgrounds into account, we obtained upper limits on the 
annihilation cross section for various models of contraction. We showed that the limits are 
improved, compared to those for the NFW or Einasto profiles with no contraction, by a factor 
of 2-6 for a wide range of parameters of the contraction model. For the canonical contraction 
model, the current limits are (o~v) < (2-3) x 10~ 25 cm 3 s" 1 for dark matter with mass around 
10 GeV, which is roughly one order of magnitude larger than the canonical cross section nec- 
essary to explain the thermal relic density. We argue that this effect is more robust than the 
subhalo boost that has been discussed in the literature, and indeed, more important unless 
the minimum subhalo mass is considerably smaller than the mass of the Sun. 
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A Analysis of Fermi-LAT data 
A.l Photon data processing 

We perform the analysis using the Fermi Science Tools (v9r23pl), as well as data of gamma- 
ray photons collected by Fermi-LAT for the mission elapsed time (MET) between 239557417 s 
and 329159098 s (2.8 years; from August 4, 2008 to June 7, 2011), contained in a region of 
interest (ROI), i.e., a circle of 10°-radius around a target cluster. First, we chose DIFFUSE and 
DATACLEAN class photons whose energies are between 1 GeV and 100 GeV by using gtselect, 
such that contamination of charged particles gets reduced to minimum. Then with gtmktime 
we further cut the data according to the spacecraft configuration, using "DATA_QUAL==1 && 
LAT_C0NFIG==1 && ABS (R0CK_ANGLE) <52" filter as well as ROI-based zenith angle cut. 

A. 2 Maximum likelihood analysis 

Using these filtered data around the cluster, we perform the binned likelihood analysis, by 
comparing the photon data with the models of sources (both point-like and extended), as 
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Figure 11. Counts map (left) and model map (right) for 1-100 GeV photons around the Fornax 
cluster. The colors show source counts per pixel and are logarithmically spaced. Ten 2FGL point 
sources are found in ROI (10° radii) and 28 total are found in the source region (15° radii), in addition 
to the both Galactic and extragalactic diffuse backgrounds. Note that the flux of the source at far 
left in the model back is saturated. 



well as the Galactic foreground and extragalactic background. The sources are taken from 
2-year source catalog (2FGL) [62], if they are located within 15° radii (source region) from 
the cluster. With gtbin, we divide the cluster map into spatial 100 x 100 pixels, whose size is 
0.2° x 0.2°, and also into 20 energy bins equally spaced logarithmically between 1 GeV and 
100 GeV. The left panel of Fig. 11 shows such a map (for Fornax) of photon counts per pixel, 
but integrated over the entire energy range, 1-100 GeV. Based on this binning, with gtltcube 
followed by gtexpcube2, we produce exposure maps corresponding to each cluster, and then 
with gtsremaps, maps of the sources and the diffuse backgrounds in ROI. 

We use gtlike to perform the maximum likelihood analysis, first with DRMNFB then fol- 
lowed by NEWMINUIT optimizers. This way, we constrain amplitudes of both the Galactic 
foreground and the extragalactic background, and also for most of the sources, the flux and 
the power-law index of the spectra if they are within ROI. If the sources are out of ROI but 
still within source region, then the fixed values from 2FGL catalog for these parameters are 
used. Then finally, by combining results of likelihood analysis, we create the model maps 
with gtmodel, to be compared with the photon count maps. The right panel of Fig. 11 shows 
the best-fit model map of the data shown in the left panel. 

A. 3 Setting upper limits on cluster emission 

To make templates of cluster emission, we convolve cluster brightness profiles with the P6_V11 
instrumental response functions (IRFs) toward the direction of ROI. Since the size of the 
cluster (i.e., r s /d,A) and the point spread function (PSF) is much smaller than ROI, we here 
use smaller regions (~ 2° x 2°) with finer pixel sizes (0.02° x 0.02°) around clusters for 
obtaining the upper limits. For this new pixelization, we again run gtbin to make photon 
count maps such as the one shown in Fig. 12, and gtmodel to make model maps of the cluster 
as well as those of the backgrounds and the other sources. We show the model maps for the 
Fornax cluster after convolving with the IRFs in Fig. 13 for both smooth NFW host-halo 
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Figure 12. The same as the left panel of Fig. 11, but with finer pixclization. White dots represent 
single photons received. 




Figure 13. The gamma-ray model maps of Fornax for the NFW density profile (left) and host plus 
subhalo model (right). The colors represent expected photon counts per 0.02° x 0.02° pixel (in linear 
scaling), where bb annihilation channel with m x = 100 GeV and (av) = 10~ 26 cm 3 s _1 is assumed, and 
photon energies are between 7.94 GeV and 10 GeV. No backgrounds and other sources are included 
in these maps. 



case and the case with large boost due to subhalos (with M m ; n = 10 _6 Mq). One can see 
that the latter is much more extended than the former. Note, however, that the backgrounds 
and other sources are not included in these maps. Although we do not show the maps for 
the contracted cluster model, they look similar to the map for NFW. 

We obtain upper limits on cluster emission by comparing the models with the data. 
According to the Bayesian statistics (e.g., [63]), the posterior probability distribution function 
for some theoretical parameter 9 (in this case, the annihilation cross section, 9 = (crv)) given 
data d is 

P(6\d) oc P(6)L(d\e), (A.l) 
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Figure 14. Model map of the Fornax cluster combined with the diffuse backgrounds as well as point 
sources, integrated for 1-100 GcV. This map is to be compared with the counts map of Fig. 12; 
the color scaling is set such that white represents single photon received per pixel, while black docs 
less than 0.007 photon. The cluster dark matter component is from the smooth NFW halo, with bb 
annihilation channel, m x = 100 GeV, and (av) = 5.5 x 10 _24 cm 3 s _1 that is the 95% credible upper 
limit from the current data. 



where P{8) is the prior distribution, for which we adopt the uniform and improper prior, 
P(8) = const., for positive 8 and zero otherwise, and L(d\8) is the Poisson likelihood function. 
The latter is specifically given as 

*■ - n < a - 2 > 

i=i 

where subscript i represents pixel number, /ij and are theoretical model and the data 
counts in the pixel i, respectively. We assume that the model is divided into the signal (sj) 
and background (bi) contributions, /ij = Sj + where the former depends on the theoretical 
parameter 8 and the latter includes both the diffuse backgrounds and 2FGL sources. We 
here fix all the parameters for bi to the values obtained with gtlike in the previous subsection, 
and then obtain an upper limit on 8 by solving 

l-a= d8 P[6\d), (A.3) 
Jo 

for 8n m . In particular, we are interested in upper limits corresponding to 95% credible 
interval, and for that we adopt a = 0.05 above. As an example, in Fig. 14, we show what 
the Fornax cluster map would appear when it is just as bright as currently allowed and is 
dominated by the smooth NFW component (with no contraction). The photon counts per 
0.02° x 0.02° pixel could be at most 0.017 at the center of the map. 

When we obtain upper limits with stacked cluster maps, we generalize Eq. (A. 2) to 
include contributions from different clusters as 
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where the subscript cl represents clusters. Here we remove some clusters from the stacking 
analysis, if they are brighter than 3a level compared with the backgrounds for any dark 
matter masses and the annihilation channels. For this purpose, we also use the posterior 
probability distribution function introduced above, P(9\d), and by fixing the parameters for 
the background models 6j. As the result, we have 38 clusters used in the stacking analysis. 
A more careful treatment for the gamma-ray detection from the galaxy clusters should be 
made by varying background parameters together with 9. However, this is beyond the scope 
of the present study and will be revisited in future publication. 
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